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ABSTRACT 


In  this  thesis,  we  study  the  feasibility  of  improving  aluminum-carbon  repulsive 
potentials  for  use  in  density-functional  tight  binding  (DFTB)  simulations  of  low-valence 
aluminum  metalloid  clusters.  These  systems  are  under  consideration  for  use  as  novel 
fuels  with  rapid  metal  combustion  kinetics,  and  contain  an  unusual  mix  of  low-valence 
metal/metal  bonds  as  well  as  organometallic  components.  We  show  that  current  DFTB 
parametrizations  of  the  repulsive  potential  for  Al/C  interactions  do  not  provide  an 
adequate  treatment  of  the  bonding  in  these  clusters.  We  performed  a  re-parametrization 
of  the  Al-C  repulsive  potential  via  comparison  to  high-level  density  functional  theory 
(DFT)  results  that  are  known  to  give  accurate  thermochemistry  for  these  clusters.  We 
found  that  the  reparametrized  system  solves  the  most  egregious  issues,  particularly  those 
associated  with  an  unphysical  distortion  of  the  q5  Al/cyclopentadienyl  bond.  DFTB 
molecular  dynamics  simulations  of  the  oxidation  of  ALjCp  4  show  reasonable  comparison 
with  a  DFT -based  Car-Parrinello  method,  including  correct  prediction  of  hydride 
transfers  from  Cp*  to  the  metal  centers  during  the  reaction. 
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I.  BACKGROUND 


A.  ROCKET  PROPELLANTS 

Solid  rocket  motors  are  widely  used  as  missile  propulsion  systems  in  the 
Department  of  Defense  (DOD)  and  space  communities.  The  typical  propellant  grain 
within  many  of  these  motors  consists  of  an  aluminum  fuel  along  with  oxidizer  and  binder 
components.  The  aluminum  particles  used  are  typically  micrometer  scale,  though  there 
has  been  considerable  recent  interest  in  aluminum  nanoparticles  to  increase  the  burning 
rate  as  well  as  lowering  the  ignition  threshold  [1].  Within  the  DOD,  rockets  are  used  in 
both  very  small  systems  such  as  the  Mk  66  solid  rocket  motor  that  propels  the  Hydra  70 
2.75  inch  rocket  and  very  large  systems  such  as  the  D-5  rocket  system  for  the  Trident  II 
missile.  As  such,  improvements  to  rocket  motor  performance,  particularly  in  areas  of 
linear  bum  rate  or  energy  density  that  may  lead  to  novel  rocket  motor  designs,  present  an 
attractive  opportunity.  Increasing  the  rate  of  aluminum  combustion,  however,  is 
challenging.  The  process  is  limited  by  the  presence  of  a  native  oxide  layer  on  the  surface 
of  the  particles,  as  well  as  the  slow  mass  transfer  that  occurs  during  its  diffusion-limited 
burning.  Recent  research  efforts  have  been  directed  at  developing  A1  nanoparticles  from  a 
different  approach,  using  solution  chemistry  methods  to  grow  small  metal  clusters  with  a 
single  monolayer  of  organic  ligand  on  their  surface.  These  may  offer  the  potential  for 
greatly  increased  combustion  kinetics,  well  beyond  that  of  typical  metal  fuels. 

B.  ALUMINUM  CY CLOPENTADIENYL  CLUSTERS 

Low-valence  A1  clusters  with  a  surface  organic  ligand  layer  have  been  known 
some  time,  discovered  in  the  pioneering  work  of  Schnockel  and  co-workers  in  recent 
decades  [2].  These  clusters  are  formed  in  a  co-condensation  reactor  which  produces 
monovalent  aluminum  halide  (A1C1  or  AlBr)  liquids  as  the  starting  material.  These 
liquids  are  then  used  for  solution  growth  of  larger  clusters,  which  includes  a  ligand 
exchange  process  with,  most  commonly,  cyclopentadienyl  (CsMes  or  Cp)  or 
hexamethyldisilazane  (funds).  During  the  growth  process,  an  unknown  kinetic 
mechanism  traps  the  system  and  prevents  a  runaway  reaction  which  would  form  zero- 
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valent  bulk  metal.  The  final  result  is  an  onion  structure  with  low-valent  metal  in  the 
interior  and  organic  ligands  on  the  outside.  These  systems  are  often  called  “metalloid,” 
since  they  contain  a  mix  of  metal-metal  bonds  as  well  as  metal/organic  interactions.  One 
of  the  most  notable  cluster  in  this  family  is  composed  of  a  50  aluminum  atoms  with  12 
Cp  ligands  strongly  bound  to  the  exterior  [2].  However,  further  refinement  and  study  of 
these  clusters  is  required  to  determine  the  optimal  balance  of  reactivity,  stability,  and 
combustion  properties  in  order  to  obtain  a  configuration  suitable  for  safe  and  effective 
use  as  a  rocket  motor  propellant  or  other  energetic  applications  [1].  The  most  efficient 
approach  to  determine  this  optimal  balance  would  be  with  the  aid  of  computational 
studies  to  guide  efforts  of  experimentalists  to  reproduce  the  desired  configurations  in 
laboratory  conditions.  Unfortunately,  further  computational  study  of  such  systems  is 
complicated  by  the  very  large  computing  resource  requirements  for  larger,  more  complex 
structures  when  using  the  current  computational  method  of  choice,  density  functional 
theory  (DFT). 

C.  DENSITY  FUNCTIONAL  TIGHT  BINDING 

Density  functional  tight  binding  (DFTB),  shows  promise  as  an  intermediate 
computational  method  that  could  yield  approximate  results  practically  close  to  the 
accuracy  of  DFT  with  significant  reduction  to  computational  resource  requirements. 
Improvements  in  the  accuracy  of  available  density  functionals,  as  well  as  the  increased 
availability  of  computational  resources,  have  made  DFT  the  preferred  computational 
method  for  electronic  structure  calculation  for  most  applications  [3].  It  provides  a  good 
balance  between  the  competing  demands  of  chemical  accuracy,  moderate  computational 
speed,  and  transferability.  However,  for  the  study  of  aluminum  nanoparticles  intended  for 
use  in  rocket  propellants,  DFT  methods  prove  to  be  quite  taxing  for  even  very  large 
supercomputers.  More  in-depth  analysis  of  systems  of  interest  requires  study  of  large 
systems,  upwards  of  500  atoms,  evolving  in  time  during  chemical  decomposition, 
oxidation,  and  similar  processes.  Many  implementations  of  DFT  do  not  scale  well  with 
the  addition  of  more  processing  units  (CPUs),  as  doubling  CPUs  assigned  may  only  yield 
a  less  than  10%  improvement  in  calculation  time  [4].  Therefore,  the  use  of  more 

approximate  methods  is  warranted  to  help  speed  research  efforts. 
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Though  derived  from  DFT,  DFTB  is  still  a  tight  binding  (TB)  method  with  the 
fundamental  assumption  of  tightly  bound  electrons.  Therefore,  it  is  ideally  suited  for 
covalent  bound  systems  such  as  hydrocarbons  [3].  Despite  the  known  limitation  of  this 
assumption,  it  has  also  shown  promising  results  in  treating  some  systems  with  ligand 
protected  metallic  clusters  (such  as  thiol-protected  gold  clusters),  with  the  benefit  of 
being  several  orders  of  magnitude  more  efficient  than  full  Kohn-Sham  DFT  computation 
[5].  However,  DFTB  requires  specific  parametrizations  to  be  performed  to  treat  each  type 
of  interaction  in  a  system,  and  the  nature  of  this  parametrization  will  greatly  affect  the 
quality  and  applicability  of  the  output.  As  discussed  below,  we  found  that  existing  DFTB 
parametrizations  in  the  literature  are  insufficient  for  study  of  the  aluminum- 
cyclopentadienyl  metalloid  cluster  systems.  The  purpose  of  this  study  is  to  develop  a  new 
repulsive  potential  parameterization  that  allows  us  to  perform  accurate  DFTB 
calculations  for  ligated  aluminum  nanoparticle  systems. 

There  are  several  DFTB  software  packages  available,  but  the  one  chosen  for  this 
research  was  DFTB+  [6].  Several  parametrization  sets,  in  the  form  of  Slater  Koster  files, 
are  also  provided  by  the  developers  free  for  academic  use.  The  parameter  set  chosen  was 
the  matsci-0-3  parameter  set,  primarily  selected  due  to  its  inclusion  of  interactions 
involving  A1  and  C,  which  appears  in  every  system  of  interest  to  this  study  [7].  The 
software’s  core  development  was  led  by  Balint  Aradi  of  the  University  of  Bremen  and 
Ben  Hourahine  of  the  University  of  Strathclyde. 
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II.  PRINCIPLES  OF  DENSITY  FUNCTIONAL  TIGHT  BINDING 


This  chapter  will  present  a  brief  overview  of  the  theory  of  DFTB.  To  maintain 
consistency  with  much  of  the  literature  that  provides  a  more  in-depth  study  of  DFTB, 
Hartree  atomic  units  will  be  used  throughout.  In  particular,  this  overview  will  draw 
heavily  from  [3],  and  all  equations  in  this  chapter  are  sourced  exclusively  from  [3]  to 
maintain  consistency  of  terms  used  throughout. 

A.  DENSITY  FUNCTIONAL  THEORY 

As  density  functional  theory  is  the  starting  point  for  DFTB,  we  must  have  at  least 
some  understanding  of  the  final  expressions  derived  in  DFT.  For  the  many  body  problem 
of  a  system  consisting  of  nuclei  and  electrons,  the  total  energy  can  be  written  as  [3]: 

E  =  T  +  Eext+Eee+En  (1) 

T  is  the  kinetic  energy,  Eext  is  the  energy  of  external  interactions,  including 
electron-ion  interactions,  Eee  is  the  electron-electron  interaction  energy,  and  Eu  is  the  ion- 
ion  interaction  energy.  In  Kohn-Sham  DFT,  with  a  system  of  non-interacting  electrons, 
the  total  energy  can  be  written  as  [3]: 

E[n(r)\  =  Ts+Eext+EH+Exc+En  (2) 

Ts  is  the  non-interacting  kinetic  energy,  Eu  is  the  Hartree  energy,  and  Exc  is  the 
exchange  correlation  (xc)  energy.  This  is  more  explicitly  written  in  Equation  3  below  [3]: 

£["]  =  Z  /.<n|(-^!+Jv„(r)n(r))|^)  +  £„[n]  +  £„  (3) 

a  ^ 

Here  fa  e  [0, 2  ]  is  the  occupation  of  a  single-particle  state  vpa  with  energy  sa.  We 
next  begin  approximating  by  considering  a  system  with  density  no{ r)  as  if  atoms  in  the 

system  were  free  and  neutral.  The  density  is  not  the  true  minimizing  density  that  gives  us 
our  ground  state  energy,  but  we  assume  it  is  close,  deviating  only  by  some  small  8n{ r) . 

A  second  order  expansion  of  Equation  3  at  na(r)  with  fluctuation  d'n(r)  provides  us  with 
Equation  4  [3]. 
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(4) 


£[*>]  =  £  f, (t"„  I~ +K„  +v„[«j+ v„[«j |r. 

a  ^ 


irrA>J,  i 


-)SnSn ' 


-JT(  .  . 

2^  SnSn'  |r-r'| 

i|v„K](r)n,(r)  +  £>„]  +  £„  -  Jv^Kl(rK(r) 


The  terms  in  Equation  4  are  normally  subcategorized  into  energy  terms.  Equation 
5  comes  from  the  first  line  in  Equation  4,  and  is  the  band  structure  energy  [3].  Equation  6 
comes  from  the  second  line,  and  is  the  energy  from  charge  fluctuations,  mainly  coming 
from  Coulomb  interactions,  but  also  containing  some  xc  contributions  [3].  Equation  7  is 
the  third  line,  and  is  called  the  repulsive  energy,  because  of  the  dominance  of  the  ion-ion 
repulsion  term,  but  also  captures  other  xc  effects  [3]. 


EBS  =Z  fa{¥a\H\-no\\Ya) 
a 


(5) 


r-r 


(6) 


^=-|jv,Kl(rK(r)  +  £„K]  +  £i,-Jv„K](rK(r)  (7) 

Equations  5-7  allow  us  to  collapse  Equation  4  into  the  more  convenient  form 
below  in  Equation  8  [3]. 

E[Sn\  =  Ebs  [Sn\  +  Ecoul  [fin]  +  Erep  (8) 

We  recall  that  atom  energy  can  be  expressed  as  a  function  of  charge  fluctuation 
Aq,  electronegativity,  and  the  Hubbard  U  parameter.  We  make  an  assumption  of  the 
Coulomb  energy  of  two  spherically  symmetric  Gaussian  charge  distributions  to  yield 
Equation  9  for  Ecoul  expressed  as  a  function  of  charge  fluctuation  [3]. 


Ecoui  =^2>z/  (Ru  )MMj 

^  ij 


(9) 
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We  then  apply  the  tight  binding  formalism  to  Equation  8,  which  are  the 
approximations  of  DFTB.  In  the  EBS  term,  we  assume  tightly  bound  valence  electrons, 

and  use  a  minimal  basis  set  with  only  one  radial  function  for  each  angular  momentum 
state.  The  Hamiltonian  matrix  elements  are  the  principal  parameters  for  this  term.  In  the 
Ecoul  term,  the  charge  fluctuations  are  approximated  with  a  Mulliken  population  analysis 
and  can  be  parametrized  by  changes  to  the  electronegativity  and  the  Hubbard  parameter 
U.  The  E  remains  the  repulsive  energy,  and  is  parametrized  for  atomic  pairs. 
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III.  FITTING  THE  REPULSIVE  POTENTIAL 


In  theory,  fitting  the  repulsive  potential  Vrep  can  be  performed  systematically  [3]. 
Vrep  is  considered  to  be  analogous  to  Exc  in  DFT,  as  it  packages  much  of  the  difficult 
physics  into  a  single  term.  However,  unlike  DFT,  DFTB  parametrization  requires 
characterizing  a  new  repulsion  for  each  pair  of  atoms  in  a  system,  the  fitting  process  can 
become  very  labor  intensive.  Additionally,  individual  parametrizations  often  require 
laborious  adjustments  that  may  not  be  fully  rigorous  with  theory.  Therefore,  for  this 
study,  we  chose  to  make  use  of  an  existing  parametrization  set,  and  limited  efforts  at 
improvement  only  for  those  parametrizations  where  we  observed  known  deficiencies. 
More  specifically,  we  attempted  to  improve  only  the  parametrization  of  the  Al-C 
interaction.  The  existing  parametrizations  that  are  used  to  treat  other  interactions  are 
contained  in  the  matsci-0-3  family  of  parametrizations  in  the  DFTB+  software  package 
[7]. 

A.  CHARACTERIZING  REPULSIVE  POTENTIAL 

The  equations  supporting  parametrization  of  the  repulsive  potential  are  as  follows 
and  are  taken  from  [3]: 

(*)  =  edft  (. R )  -  [Ebs  (R) + Ecoul  (*)]  (10) 

N  •  Vrep  (R)  =  EDFr(R)~  [Ebs  (R)  +  Ecou,  (*)]  (11) 

Equation  10  is  used  for  a  system  with  a  single  element  pair.  Equation  11  is  used 
for  a  symmetric  structure  where  there  are  N  bonds.  In  practice,  a  single  system  is 
insufficient  to  provide  a  repulsive  potential  that  can  be  robustly  applied  to  a  variety  of 
different  systems  with  different  bonding  types.  In  the  existing  matsci-0-3  parametrization 
set,  the  Al-C  interaction  was  created  for  the  investigation  of  the  interaction  of 
ethylphosphonic  acid  (C2H5PO(OH)2)  with  aluminum  oxide  surfaces  and  was  fit  using  a 
CH3-AIH2  system  [8].  As  this  system  consists  of  a  different  bonding  type  from  the 
delocalized  bonding  present  in  ligand  protected  metallic  clusters,  it  would  have  been  a 
reasonable  assumption,  even  without  test  calculations,  that  a  re-parametrization  would  be 
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necessary.  Unsurprisingly,  we  observed  in  test  calculations  that  re-parametrization  is 
indeed  required,  as  the  test  structure  (Figure  1)  shows  a  significant  and  obvious  distortion 
in  the  geometry  of  the  AlCp  monomer  that  is  a  building  block  to  our  systems  of  interest. 


Figure  1 .  AlCp  Monomer 

Therefore,  we  start  with  relatively  simple  structures  of  interest,  where  DFT  can 
still  quickly  provide  us  with  geometries  R,  energies  E\wi(R),  and  forces  Fdft  (which 
would  be  zero  for  an  optimized  structure).  One  approach  would  be  to  fit  Vrep  to  minimize 
the  energy  and  force  differences  and  l£m  1  -  £,DFTB|  and  IFdft  -  F,dftbI.  However,  as 
advised  by  [3],  in  practice,  we  only  minimize  the  force  difference.  The  reason  for  this 
recommendation  is  that  it  provides  us  the  benefit  of  working  with  absolute  terms  as, 
unlike  energies,  forces  are  absolute  terms.  Additional  justification  is  provided  in  that 
most  of  the  energy  contribution  in  DFTB  comes  from  the  band  structure  term.  This  term 
is  adjusted  by  modifying  the  confinement  potential  and  the  Hubbard  parameter.  If  we  find 
that  we  need  to  make  adjustments  to  minimize  the  energy  difference  at  large  distances, 

we  cannot  fix  the  errors  in  energy  by  adjusting  the  short-range  repulsion  energy  alone. 
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Additionally,  repulsion  energy  needs  to  be  monotonic  and  smooth.  If  we  adjust  repulsion 
too  abruptly,  the  parametrization’s  forces  will  suffer  from  the  erratic  repulsion.  For  the 
purposes  of  this  study,  we  did  not  make  any  adjustments  to  the  confinement  potential  or 
Hubbard  parameter.  Therefore,  the  fit  is  only  to  the  forces,  or  the  derivative  of  the 
repulsive  potential,  yielding  [3]: 

1/  _  EDft  (RAB  )  _  I  RrS  (RaB  )  +  R<:oul  ( RAB  )]  /  1  O  \ 

^  ~  N  K  ’ 

B.  FITTING  THE  REPULSION 

The  systems  chosen  as  the  base  of  our  parametrization  effort  were  the  simple 
CH3-AIH2  system  used  in  the  original  matsci-0-3  parametrization,  a  trimethylaluminum 
(A1(CH3)3)  system,  and  an  aluminum  cyclopentadienyl  (AlCp)  monomer.  We  can  rapidly 
generate  V  rep  curves  for  these  systems  by  using  DFT  and  DFTB  to  obtain  EDFTt  EBs,  and 
Ecol,i  by  varying  the  distance  RAb  to  obtain  a  range  of  energy  values  as  a  function  of  RAB. 

1.  Collecting  Data 

DFT  calculations  to  obtain  EDFT  for  these  relatively  simple  systems  were 
performed  by  optimizing  the  geometries  using  the  M06-2X  functional  and  6-3 1G*  basis 
set  in  the  Gaussian  09  software  package.  This  combination  of  functional  and  basis  set  had 
been  previously  determined  by  our  research  group  to  give  accurate  bonding  energies  of 
aluminum  metalloid  clusters  (typically  within  5  kJ/mol  of  experiment). 

The  energy  calculations  to  obtain  EBS  and  Ecou/  were  performed  using  DFTB+  on 
the  same  geometries  with  varying  RAB  as  those  used  in  the  DFT  calculations.  All  DFTB+ 
calculations  were  performed  using  Version  1.2  of  the  DFTB+  software  package,  using 
Parser  Version  4.  Repulsive  interaction  was  characterized  using  the  polynomial  option. 
The  Hamiltonian  was  calculated  using  SCC  (self-consistent  charge)  calculations,  with  an 

o 

SCC  tolerance  of  1.0E-6  e,  and  the  maximum  force  component  was  0.02  eV/  A. 
Derivatives  EDFT  ,  EBS  ,  and  Ecoul  were  obtained  with  numerical  differentiation. 
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2.  Fitting  Equations 

The  DFTB+  software  package  allows  for  a  polynomial  characterization  of  the 
repulsive  potential  as  shown  in  Equation  13  [9]: 

v„=2>,-  <13> 

i=2 

From  this  form  of  V  ,  we  can  easily  obtain  an  expression  for  the  derivative  of 
the  repulsive  potential,  with  which  we  can  then  fit  our  data. 

v'reP  =-Yjici(rcu,-ry~1  (14) 

1=2 


(15) 


(16) 


As  the  force  is  the  negative  of  the  derivative  of  potential  energy,  the  coefficients 
as  determined  from  a  fit  to  the  force  can  then  be  directly  used  in  the  polynomial  equation 
for  the  repulsive  potential.  An  example  of  a  fitted  force  curve  from  previous  literature  on 
the  boron-nitrogen  interaction  is  provided  in  Figure  2  [10]. 
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Figure  2.  Example  Force-Based  Fitting  of  the  Boron-Nitrogen  Repulsion 

Parametrization 

It  will  be  clear  to  the  most  casual  observation  that  this  fitting  is  not  in  particularly 
good  agreement  with  the  plotted  points.  We  will  delve  into  the  reasons  for  this  in  our 
discussion  of  the  fitting  procedure. 

3.  Fitting  Procedure 

An  initial  test  fitting  was  made  solely  to  AlCp  using  a  least  squares  fitting 
algorithm  to  the  form  of  the  polynomial  in  Equation  16.  For  the  value  of  rcut,  the  original 
parametrization ’s  value  of  6.81  Bohr  was  left  unchanged.  Unfortunately,  initial  attempts 
using  just  the  output  of  the  algorithm  were  not  wholly  adequate.  While  fitting  alone  was 


13 


able  to  provide  a  result  that  corrected  the  ring  slip  deficiency  caused  by  the  stock  matsci- 
0-3  parametrization  in  the  AlCp  monomer,  it  provided  Vrep  curves  that  were  not  well 
aligned  with  physical  conventions.  Therefore,  small  manual  adjustment  of  the  polynomial 
coefficients  was  conducted  with  performance  testing  on  a  family  of  AlCp  based  systems 
with  a  focus  on  improving  the  geometry  results.  Some  compromise  was  required  between 
the  accuracy  of  the  calculation  and  smoothness  and  monotonicity  of  the  function  in  order 
to  avoid  compromising  transferability.  Specifically,  in  addition  to  the  monomer, 
adjustments  were  tested  against  both  simple  covalent  bonding  systems  and  tetramer 
geometries.  The  tetramers  proved  to  be  most  sensitive  to  adjustments  to  small  changes  to 
the  Vrep  of  the  system,  with  even  small  changes  leading  to  distortions  in  the  Al-Al 
bonding  of  those  clusters.  During  this  adjustment,  it  was  observed  that  changes  in  the 
region  of  interest  (between  4-5  Bohr)  provided  the  most  significant  impact  to  the 
accuracy  of  output  geometries  for  this  family  of  systems. 


Figure  3  shows  the  fit  curve  for  the  force,  which  is  plotted  along  with  force  curves 
obtained  for  the  three  systems  discussed  in  Section  1,  for  which  we  obtained  DFT  data. 
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Visually,  the  force  curves  are  not  in  particularly  good  agreement  with  any  of  the  systems. 
The  fitting  algorithm  alone  did  not  produce  satisfactory  results  when  all  test  geometries 
were  run,  therefore  manual  adjustments  were  necessary.  Manual  adjustments  were 
performed  on  Vrep,  not  on  the  force  curves,  therefore  the  fit  to  the  force  suffered. 
However,  as  seen  previously  in  Figure  2,  this  kind  of  result  is  not  particularly  unusual  in 
DFTB  potentials.  In  the  parametrization  in  Figure  2,  the  author  chose  to  make 
adjustments  to  a  Vrep  curve  from  a  strict  fit  to  the  data  in  order  to  match  the  results  of 
band  structure  calculations  for  his  systems. 

For  our  parametrization,  transferability  to  our  systems  of  interest  was  the 
paramount  concern.  A  spline  is  often  used  in  the  fitting  of  other  DFTB  parametrizations, 
which  would  provide  the  benefit  of  being  able  to  more  finely  adjust  fitting  for 
experimentation  with  parametrization  performance  to  isolate  changes  to  more  specific 
ranges  of  separations.  However,  we  chose  the  polynomial  fit  due  to  its  easier 
implementation  of  re-parametrization  efforts  in  DFTB+.  Table  1  contains  the  final 
parametrization ’s  coefficients: 

Table  1.  Vrep  Polynomial  Coefficients 


Coefficient 

Value 

c2 

0.0352 

c3 

-0.067 

c4 

0.0585 

c5 

-0.0273 

c6 

0.0069 

c7 

-0.0010 

c8 

0.0001 

15 


As  the  coefficients  obtained  by  fitting  the  force  are  identical  to  the  coefficients  in 
our  repulsive  potential,  we  can  then  plug  these  coefficients  in  directly  to  Equation  13, 
which  provides  us  with  our  repulsive  potential  shown  in  Figure  4. 


Figure  4.  Repulsive  Potential  Vrep  from  This  Work  and  the  Matsci-0-3  Set. 

Notable  features  in  Figure  4  are  that  the  new  repulsive  potential  is  significantly 
more  repulsive  for  distances  that  are  shorter  than  those  typically  encountered  in  the 
systems  we  are  interested  in,  which  is  approximately  where  the  two  curves  intersect,  and 
roughly  coincides  with  the  Al-C  separation  found  in  most  DFT  equilibrium  geometry 
calculations  for  the  AlCp  monomer.  As  we  increase  to  distances  greater  than  that 
equilibrium  separation,  the  repulsive  potential  becomes  significantly  lower. 
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IV.  TESTING  THE  NEW  REPULSIVE  POTENTIAL 


We  begin  testing  our  new  repulsive  potential  with  an  examination  of  the  output 
geometries.  DFT  calculations  for  these  systems  were  run  using  the  Gaussian  09  software 
package,  using  the  6-3 1G*  basis  set  and  M06-2X  functional.  The  only  exception  was  the 
Al5oCp*i2  cluster,  which  was  optimized  using  the  B3LYP  functional.  DFTB  calculations 
were  run  with  DFTB+  Version  1.2,  using  both  the  stock  matsci-0-3  parameter  set  and  a 
modified  matsci-0-3  parameter  set  using  the  reparametrized  Al-C  interaction. 

A.  GEOMETRY  OPTIMIZATION 

In  addition  to  a  qualitative  examination  of  the  geometries  obtained  through 
DFTB,  we  also  performed  measurements  of  the  bond  lengths  to  perform  a  quantitative 
examination  of  our  new  repulsive  potential’s  performance.  We  also  measure  the  q 
parameter,  which  quantifies  how  well  the  bonding  of  the  A1  atom  is  centered  on  the  Cp 
ring.  The  q  parameter  is  also  considered  to  be  represent  how  many  of  the  ring  C  atoms 
are  bonded  to  the  central  A1  atom  [11].  An  q5  value  would  correspond  to  equal  bonding 
between  the  A1  atom  and  each  of  the  five  C  ring  atoms,  while  an  q1  value  would 
correspond  to  bonding  predominantly  between  the  A1  atom  and  just  one  of  the  five  C  ring 
atoms. 


1.  Monomers 

The  first  system  is  the  same  that  was  used  to  demonstrate  the  deficiencies  in  the 
current  Al-C  parametrization  in  the  DFTB+  matsci-0-3  parameter  set.  The  AlCp 
monomer  is  the  most  basic  building  block  for  the  ligand  protected  aluminum  clusters  we 
are  interested  in.  Figure  4  shows  the  geometries,  with  Figures  5a  and  5d  depicting  the 
structure  optimized  with  DFT,  Figure  5b  and  5e  depicting  the  structure  optimized  with 
the  original  DFTB  parameter  set,  and  Figure  5c  and  5f  depicting  the  structure  optimized 
with  the  reparametrized  DFTB  parameter  set. 
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Boxes  a)  and  d)  depict  DFT  optimization.  Boxes  b)  and  e)  depict  stock  DFTB+ 
optimization.  Boxes  c)  and  f)  depict  reparametrized  DFTB+  optimization.  Note  the 
shifting  of  the  A1  atom  and  distortion  in  the  Cp  ring  present  in  stock  DFTB+  optimization 
that  is  corrected  in  the  reparametrized  DFTB  optimization. 

Figure  5.  AlCp  Monomer — Side  and  Head-on  Views 


A  qualitative  visual  inspection  of  the  geometries  reveals  that  the  optimized 
geometry  produced  by  a  reparametrized  DFTB  calculation  provides  a  reasonably  close 
approximation  to  the  DFT  geometry.  The  visually  apparent  distortion  due  to  ring  slip  in 
the  bonding  between  the  A1  atom  and  the  Carbon  ring  in  the  AlCp  molecule  is  no  longer 
present,  as  the  A1  atom  is  properly  centered  in  the  q5  position. 

Figure  6  depicts  an  AlCp*  monomer,  which  is  a  methylated  AlCp,  with  a  CH3 
methyl  group  replacing  each  H  atom  attached  to  each  C  atom  on  the  AlCp  carbon  ring.  It 
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is  of  greater  interest  as  a  component  of  more  complex  ligand  protected  aluminum  clusters 
that  hold  promise  in  energetics  applications.  Figure  6a  and  6d  present  the  DFT  optimized 
structure,  Figures  6b  and  6e  present  the  results  of  the  original  DFTB  parameter  set,  and 
Figures  6c  and  6f  are  from  the  reparametrized  DFTB  parameter  set. 


AlCp* 


Matsci-0-3  DFTB+ 


Reparametrized  DFTB 


AlCp* 


DFT 


Matsci-0-3  DFTB+ 


Reparametrized  DFTB 


Boxes  a)  and  d)  depict  DFT  optimization.  Boxes  b)  and  e)  depict  stock  DFTB+ 
optimization.  Boxes  c)  and  f)  depict  reparametrized  DFTB+  optimization. 
Reparametrized  interaction  is  again  necessary  to  correct  the  At  atom  shift  and  Cp  ring 
distortion. 

Figure  6.  AlCp*  Monomer — Side  and  Head-on  Views 


The  reparametrized  DFTB  result  is  again  reasonably  close  visually  to  the  DFT 


result,  and  offers  a  significant  improvement  over  the  original  DFTB  parametrization, 
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which,  similar  to  the  AlCp  results,  shifts  the  A1  atom  significantly  from  the  correct  r|5 
position. 


Box  a)  depicts  DFT  optimization.  Box  b)  depicts  stock  DFTB+  optimization.  Box  c) 
depicts  reparametrized  DFTB+  optimization.  Surprisingly,  the  stock  DFTB+  parameter 
set  does  not  cause  clearly  apparent  distortion  for  this  particular  monomer. 

Pr 

Figure  7.  AlCp*  Monomer — Head-on  Views 


Pr 

Figure  7  depicts  an  AlCp*  monomer,  which  removes  one  of  the  methyl  groups 
from  the  AlCp*  monomer,  and  replaces  it  with  a  propyl  group.  In  these  views,  the  propyl 
group  is  attached  to  the  top  of  the  ring.  Once  again,  the  reparametrized  DFTB  provides  a 
visually  satisfying  approximation  of  the  DFT  result.  What  is  surprising  is  the  original 
DFTB  result  also  does  not  present  any  obvious  distortion  to  the  geometry  from  the  DFT 
result.  However,  as  stated  previously  when  discussing  the  fitting  strategy,  visually 
satisfactory  results  in  monomer  optimization  did  not  necessarily  carry  over  to  the  next 
level  in  complexity  for  these  structures,  as  will  be  discussed  further  in  the  next  section. 

2.  Tetramers 

Figure  8  shows  an  AlCp  tetramer  (AI4CP4),  formed  by  combining  four  AlCp 
monomers  with  the  core  A1  atoms  in  a  tetrahedral  shape.  The  deficiency  in  the  original 
DFTB  parametrization  is  very  apparent  in  this  view,  as  the  core  A1  atoms  maintain  the 
proper  configuration,  the  Cp  rings  both  tilt  and  shift  from  their  correct  positions,  and 
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there  is  also  noticeable  distortion  in  the  Cp  ring,  as  the  H  atom  closest  to  the  A1  atom 
bonded  to  the  ring  begins  to  tilt  its  bond  out  and  away  the  ring  configuration. 


Reparametrized  DFTB 


DFT 


AI4Cp4 


Matsci-0-3  DFTB+ 


Box  a)  depicts  DFT  optimization.  Box  b)  depicts  stock  DFTB+  optimization.  Box  c) 
depicts  reparametrized  DFTB+  optimization.  This  tetramer  exhibits  a  carryover  of  the 
distortion  present  in  its  component  monomers.  The  reparametrized  DFTB  result  closely 
approximates  the  DFT  result. 

Figure  8.  AlCp  Tetramer  (AI4CP4) 


Similar  issues  are  apparent  in  Figure  9,  which  depicts  an  AlCp*  tetramer 
(Al4Cp*4).  Here  the  Cp  ring  also  tilts  and  shifts,  with  the  same  distortion  in  the  Cp  ring, 
except  in  this  case,  the  methyl  group,  which  has  replaced  the  H  closest  to  the  A1  atom, 
tilts  away  from  the  ring.  The  reparametrized  DFTB  result  significantly  improves  upon 
these  deficiencies,  and  appears  to  be  in  good  agreement  with  DFT  results. 
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Matsci-0-3  DFTB+ 


a)  AI4Cp*4 


DFT 


Reparametrized  DFTB 


Box  a)  depicts  DFT  optimization.  Box  b)  depicts  stock  DFTB+  optimization.  Box  c) 
depicts  reparametrized  DFTB+  optimization.  As  in  the  previous  image,  the  tetramer 
carries  over  the  distortion  of  its  component  monomers,  which  is  corrected  by  the 
reparametrized  DFTB  interaction. 

Figure  9.  AlCp*  Tetramer  (ALjCp*^ 


Pr 

The  final  tetramer  Al4Cp*  4  reveals  similar  issues  as  the  other  tetramer  systems. 
Significant  shift  is  evident  in  the  ligand  positioning  from  the  proper  positions  around  the 
A1  core  in  the  stock  DFTB+  optimization  depicted  in  Figure  10b.  This  occurs  despite  the 
visually  satisfying  monomer  in  Figure  7b.  Distortion  of  the  ligand  itself  is  present  as  well, 
with  the  Propyl  groups  being  tilted  out  and  away.  This  further  validates  a  need  to  extend 
some  testing  to  more  complicated  systems  before  final  selection  of  a  parametrization  of 
the  repulsive  potential.  Once  again,  the  reparametrized  potential  corrects  the  visually 
evident  deficiencies. 
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Box  a)  depicts  DFT  optimization.  Box  b)  depicts  stock  DFTB+  optimization.  Box  c) 
depicts  reparametrized  DFTB+  optimization.  Despite  the  visually  satisfactory  result  in 
the  component  monomer,  the  stock  parametrization  leads  to  the  same  category  of 
distortion  in  the  other  tetramers  in  this  family  of  systems. 

Figure  10.  AlCp*Pl  Tetramer  (ALjCp*1^) 


3.  Systems  of  Greater  Complexity 

Unfortunately,  DFTB  is  still  approximate  relative  to  DFT,  and  with  a  very  limited 
scope  of  change,  we  encounter  some  trouble  with  more  complex  systems.  The  first 
system  with  an  added  layer  of  complexity,  Al8Cp*4,  is  depicted  in  Figure  11,  and  is 
closely  related  to  the  tetramer  AUCp*4.  It  adds  4  A1  atoms  to  the  A1  core  in  the  tetramer 
Al4Cp*4. 
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a)  AI8Cp*4 

b) 

c) 

DFT 

Matsci-0-3  DFTB+ 

cm  cj^ 

Reparametrized  DFTB 

Box  a)  depicts  DFT  optimization.  Box  b)  depicts  stock  DFTB+  optimization.  Box  c) 
depicts  reparametrized  DFTB+  optimization.  The  stock  DFTB+  parametrization  exhibits 
significant  asymmetric  shifting  of  the  A1  core  in  addition  to  the  shifting  and  distortion  of 
the  Cp*  ligands.  However,  the  reparametrized  DFTB  result  causes  some  distortion  in  the 
Al-Al  core  as  well. 

Figure  11.  AlgCp*4 

Unsurprisingly,  the  stock  DFTB  parametrization  causes  similar  distortion  in  both 
the  Cp*  ligand  positioning  about  the  core,  as  well  as  within  the  ligands  themselves.  And 
again,  the  reparametrized  DFTB  result  correctly  positions  the  Cp*  ligands  about  the  core, 
without  any  apparent  distortion  in  the  ligands.  However,  we  do  see  some  differences  in 
the  Alg  core.  In  the  DFT  optimized  structure,  we  see  that  the  AF  core  is  a  single 
tetrahedron  of  4  aluminum  atoms,  with  an  inverted  tetrahedron  that  is  contained  inside 
the  larger  outer  tetrahedron.  In  our  new  DFTB  result,  the  inner  tetrahedron  expands 
outside  the  boundary  formed  by  the  outer  tetrahedron. 

A  system  of  great  interest  and  of  far  greater  complexity  is  AFoCp*i2,  depicted  in 
Figure  12.  This  system  is  the  type  of  cluster  with  potential  for  use  in  rocket  fuels  and 
other  energetics,  as  it  contains  a  sizable  internal  aluminum  nanocluster  with  protective 
ligands. 
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Al5C)Cp*i2 


DFT 


Matsci-0-3  DFTB+ 


Reparametrized  DFTB 


Box  a)  depicts  DFT  optimization.  Box  b)  depicts  stock  DFTB+  optimization.  Box  c) 
depicts  reparametrized  DFTB+  optimization.  The  stock  DFTB+  parametrization  exhibits 
significant  asymmetric  shifting  within  the  A1  core  in  addition  to  shifting  and  distortion  of 
the  Cp*  ligands.  There  is  no  obvious  distortion  in  the  A1  core  with  the  reparametrized 
DFTB  result,  however,  some  Cp*  ligands  shift  from  their  DFT  positions. 


Figure  12.  Al50Cp*i2 


The  resulting  calculations  again  demonstrate  some  of  the  limitations  of  the 
approximate  method.  While  it  does  improve  visibly  upon  the  results  of  the  stock  DFTB 
parametrization,  it  is  not  able  to  match  the  DFT  result  completely.  Specifically,  the  Cp* 
rings  in  the  new  parametrization  begin  to  demonstrate  a  slight  tilt,  though  the  positioning 
of  the  ligands  themselves  is  not  significantly  inferior.  Additionally,  the  new 
parametrization  improves  upon  the  significant  distortion  to  the  internal  A1  core  present  in 
the  stock  parametrization. 


4.  Bond  Lengths 

A  comparison  of  bond  lengths  provides  us  with  a  more  quantitative  means  of 
evaluating  the  accuracy  of  DFTB  relative  to  DFT.  The  bond  lengths  of  interest  are  the  C- 
C  distances  (in  the  Carbon  ring  in  the  Cp  component),  the  Al-C  distances  (for  more 
complex  structures,  between  the  Cp  ring  and  its  bonded  A1  atom),  the  Al-Al  distances 
(where  applicable)  and  the  C-H  distances.  We  are  also  interested  in  the  Ring  Slip 
distance,  which  is  defined  to  be  the  distance  between  the  perpendicular  projection  of  the 
bonded  A1  atom  onto  the  plane  of  its  corresponding  C5  ring  and  the  center  of  the  C5  ring. 
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This  distance  is  used  to  determine  the  p  parameter,  or  hapticity,  which  has  been 
categorized  by  [10]  as:  p5  (0  A),  p3  (0.8  A),  p2  (1.0  A),  and  p1  (>1.2  A). 

Table  2.  Bond  Lengths  (DFT) 


Cluster 

Al-C 

(A) 

Al-Al 

(A) 

C-C  Ring 
Carbons 
(A) 

C-H(A) 

A1H2CH3 

1.953 

- 

- 

1.094 

Tri- 

Methyl- 

Aluminum 

1.960 

- 

- 

1.094 

AlCp 

2.35 

- 

1.413 

1.079 

AlCp* 

2.32 

- 

1.424 

1.095 

AlCp*Pr 

2.325 

- 

1.424 

1.095 

AI4CP4 

2.357 

2.721 

1.416 

1.082 

Al4Cp*4 

2.331 

2.768 

1.423 

1.096 

Al4Cp*Pr4 

2.335 

2.783 

1.423 

1.096 

Al8Cp*4 

2.250 

2.629 

1.427 

1.095 

Table  3.  Bond  Lengths  (Reparametrized  DFTB) 


Cluster 

Al-C 

(A) 

Al-Al 

(A) 

C-C  Ring 
Carbons 
(A) 

C-H(A) 

A1H2CH3 

2.261 

- 

- 

1.105 

Tri- 

Methyl- 

Aluminum 

2.264 

- 

- 

1.106 

26 


AlCp 

2.46 

- 

1.425 

1.097 

AlCp* 

2.48 

- 

1.434 

1.112 

AlCp*Pr 

2.48 

- 

1.434 

1.112 

AI4CP4 

2.49 

2.807 

1.425 

1.097 

Al4Cp*4 

2.51 

2.806 

1.433 

1.115 

Al4Cp*Pr4 

2.51 

2.806 

1.433 

1.115 

AlsCp*4 

2.5 

2.768 

1.433 

1.113 

Table  4.  Bond  Lengths  (Stock  DFTB+  Matsci)  -  In  Progress 


Cluster 

Al-C 

(A) 

Al-Al 

(A) 

C-C  Ring 
Carbons 
(A) 

C-H(A) 

A1H2CH3 

1.970 

- 

- 

1.110 

Tri- 

Methyl- 

Aluminum 

1.980 

- 

- 

1.108 

AlCp 

3.354 

- 

1.437 

1.104 

AlCp* 

3.087 

- 

1.216 

1.112 

AlCp*Pr 

2.705 

- 

1.432 

1.113 

Al4Cp4 

3.351 

2.769 

1.033 

1.103 

Al4Cp*4 

3.058 

2.775 

1.215 

1.115 

Al4Cp*Pr4 

3.049 

2.769 

1.219 

1.114 

Al8Cp*4 

3.3072 

2.775 

1.446 

1.112 
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Table  5.  Ring  Slip  Parameter  q 


Cluster 

DFT, 

n 

parameter, 

(A) 

Reparametrized 
DFTB, 
ij  parameter, 
(A) 

Stock 

matsci-0-3 

DFTB, 

0 

parameter, 

(A) 

AlCp 

0.00017 

Ol5) 

0.0001  Ol5) 

2.825  (if) 

AICp* 

0.00022 

in5) 

0.0011  (T|5) 

0.963  (if) 

AlCp*Pr 

0.00703 

Ol5) 

0.0013  (if) 

0.0252  (if) 

Al4Cp4 

0.00307 

Ol5) 

0.0166  0i5) 

2.405  (if) 

Al4Cp*4 

0.04801 

in5) 

0.0335  (if) 

2.678  O11) 

AI4Cp*Pr4 

0.06371 

Ol5) 

0.025  (if) 

1.032  (if) 

Al8Cp*4 

0.38842 

in5) 

0.0031  (if) 

- 

All  reparametrized  DFTB  and  DFT  results  have  if  values  of  hapticity. 

The  ring  slip  for  the  stock  DFTB  matsci  parametrization  produces  unacceptably  large  shifts  in  r| 
value. 


The  quantitative  data  in  Tables  2-5  supports  the  qualitative  assessment  of  the 
geometry  visualizations.  The  stock  matsci-0-3  parametrization  demonstrates  significant 
problems  throughout  the  family  of  aluminum  cyclopentadienyl  clusters.  More 

o 

specifically,  we  see  that  the  Al-C  bond  lengths  are  significantly  longer  (almost  1  A  in  the 
worst  case)  than  our  DFT  results.  We  can  also  quantify  the  ring  slip  and  q  parameter, 
with  all  but  one  monomer  shifting  to  the  q  or  q  position.  However,  it  does  performs 
well  for  simpler  C-Al  bonds,  such  as  those  found  in  the  CH3-AIH2  (used  in  the  original 
matsci-0-3  parametrization)  and  trimethylaluminum  (A1(CH3)3)  systems. 

The  Al-C,  Al-Al,  C-C,  and  C-H  distances  of  both  the  reparametrized  DFTB 
results  and  the  DFT  results  show  good  agreement.  C-C  and  C-H  bond  lengths  for  clusters 

o 

computed  with  reparametrized  DFTB  and  DFT  are  less  than  0.01  A  apart,  Al-Al 
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o  o 

distances  are  well  within  0.1  A,  and  C-Al  distances  are  within  0.2  A.  However,  we  lose 
some  transferability  to  more  conventional  Al-C  bonds,  with  CH3-AIH1  and 

o 

trimethylaluminum  (A1(CH3)3)  bond  lengths  becoming  approximately  0.3  A  longer  than 
DFT  or  the  stock  DFTB  results.  However,  this  is  an  acceptable  penalty  for  study  of  the 
systems  we  are  interested  in,  as  these  types  of  Al-C  bonds  are  not  expected  in  the  reacted 
products. 

B.  ELECTRONIC  STRUCTURE  CALCULATIONS 

The  HOMO-LUMO  gap  is  the  energy  difference  between  the  highest  occupied 
molecular  orbital  (HOMO)  and  lowest  unoccupied  molecular  orbital  (LUMO).  We  can 
see  in  Table  6  that  the  DFTB  result  consistently  underestimates  the  HOMO-LUMO  gap 
in  comparison  to  the  DFT  result.  Interestingly,  the  amount  by  which  DFTB 
underestimates  the  HOMO-LUMO  gap  is  consistent  within  each  family  of  systems.  For 
the  monomers,  the  underestimation  is  approximately  0.8  eV,  while  the  tetramers 
underestimate  by  approximately  2.0  eV.  The  AlsCp*4,  which  is  somewhat  structurally 
similar  to  the  tetramers  has  a  difference  of  approximately  2.6  eV. 

Table  6.  HOMO-LUMO  Gap 


Cluster 

DFT 

(eV) 

DFTB 

(eV) 

AlCp 

5.869 

5.090 

AlCp* 

5.624 

4.874 

AlCp*Pr 

5.616 

4.866 

Al4Cp4 

4.832 

2.831 

Al4Cp*4 

4.626 

2.623 

Al4Cp*Pr4 

4.586 

2.663 

Al8Cp*4 

3.643 

1.031 

This  result  is  not  unreasonable,  as  the  Slater-Koster  files  for  the  matsci  parameter 
set  for  these  elements  indicate  that  it  was  originally  built  using  the  LDA  functional.  The 

HOMO-LUMO  gap  for  our  DFT  calculations  were  performed  with  the  M062X 
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functional,  and  comparison  testing  indicated  that  LDA  functionals  also  underestimated 
the  HOMO-LUMO  gap  in  comparison  for  these  systems  relative  to  the  M062X 
functional.  Recall  from  Equation  12  that  DFT  energy  is  key  to  building  the  DFTB 
parametrization.  Therefore,  we  expect  that  this  underestimation  of  the  HOMO-FUMO 
gap  is  inherited  from  the  FDA  functional  the  matsci  parameter  set  was  originally  built  on. 
However,  FDA  performs  well  for  treatment  of  ligand  interaction  with  bulk  materials  [12]. 
Therefore,  it  was  believed  that  continuing  to  base  our  parametrizations  on  the  FDA 
functional  would  provide  better  transferability  to  larger  clusters. 


Table  7 .  B  inding  Energy 


Cluster 

DFT 

(kj/mol) 

DFTB 

(kj/mol) 

Measured 

(kj/mol) 

Al4Cp4 

140 

98.8 

- 

Al4Cp*4 

159 

142 

150 

Al4Cp*Pr4 

167 

143 

160 

The  binding  energy  of  the  tetramers  is  calculated  by  taking  the  difference  between 
the  total  energy  of  the  tetramer  and  the  total  energy  of  four  of  the  component  monomers, 
and  is  presented  in  Table  7.  The  DFTB  result  shows  that  the  tetramers  are  slightly  weakly 
bound  in  comparison  to  DFT.  The  DFT  M062X  functional  was  specifically  chosen  to 
give  good  agreement  between  DFT  calculated  binding  energy  values  and  two  known 
experimental  values  [13].  For  the  two  experimentally  validated  systems,  we  have 
generally  good  agreement  between  DFT  and  DFTB.  Additionally,  the  magnitude  of 
difference  between  our  reparametrized  DFTB  results  and  DFT  results  are  significantly 
better  than  that  which  was  seen  with  some  other  DFT  functionals. 
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Table  8.  Dipole  Moment 


Cluster 

DFT 

(Debye) 

DFTB 

(Debye) 

X 

y 

Z 

total 

X 

y 

z 

total 

AlCp 

1.126 

-0.0006 

-0.0007 

1.13 

0.747 

-0.00061 

-0.00104 

0.747 

AlCp* 

-0.001 

0.0008 

-1.780 

1.78 

-0.0618 

0.120 

-0.381 

0.404 

AlCp*Pr 

-0.470 

-0.0001 

-1.791 

1.85 

0.0952 

-0.00000062 

0.364 

0.376 

Al4Cp4 

-0.0051 

0.0075 

-0.0154 

0.0179 

-0.110 

0.0389 

0.0602 

0.131 

Al4Cp*4 

0.0519 

0.0705 

-0.134 

0.160 

-0.307 

-0.0237 

-0.131 

0.334 

Al4Cp*Pr4 

-0.0289 

-0.0104 

-0.0798 

0.0855 

0.0199 

0.0311 

0.0141 

0.0396 

Al8Cp*4 

0.0389 

-0.0068 

0.0255 

0.0470 

0.00628 

0.00497 

-0.00312 

0.0086 

Dipole  moment  calculations  from  our  DFTB  results  present  some  significant 
discrepancies  in  the  monomers.  As  we  see  in  Table  8,  for  the  monomers,  DFT  is  in  good 
agreement  with  qualitative  expectations,  aligning  the  dipole  moment  to  have  a  strong 
component  along  whichever  axis  is  passing  through  the  center  of  the  Cp  ring.  Though 
this  is  still  somewhat  present  in  the  DFTB  results,  both  the  total  magnitude  and 
directionality  of  the  dipole  moment  is  not  in  good  agreement  with  the  DFT  results.  This 
becomes  particularly  problematic  for  the  AlCp*  and  AlCp*Pl  monomers.  This  is 
somewhat  consistent  with  plots  of  calculated  charge  density  for  these  monomers,  where 
we  observed  asymmetric  charge  densities  around  the  Cp  ring,  clustering  a  concentrating  a 
higher  charge  density  for  one  of  the  C  atoms  in  the  Cp  ring,  which  is  not  the  expected 
result. 

However,  the  performance  improves  for  the  tetramers  (including  AlsCp*,*),  where 

we  would  expect,  as  is  predicted  by  DFT,  that  the  total  dipole  moment  would  be  rather 
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small,  as  the  systems  generally  exhibit  a  symmetrical  configuration  about  each  core.  Here 
there  is  better  agreement  between  the  DFTB  and  DFT  results,  with  less  than  order  of 
magnitude  difference  between  the  total  dipole  moment  values. 

C.  MOLECULAR  ORBITALS  AND  ENERGY  LEVELS 

The  DFTB+  software  package  includes  a  molecular  orbital  plotting  utility  called 
waveplot.  We  used  this  utility  to  plot  the  molecular  orbitals  for  monomers  and 
investigated  its  applicability  for  tetramers  as  well. 


Box  a)  depicts  the  HOMO  for  AlCp  and  its  nonbonding  interaction  between  the  Cp  a  , 
and  the  At  sp ,  Box  b)  depicts  the  HOMO-7  for  AlCp  and  its  Cp  a  t  bonding  interaction 
with  the  At  sp,  and  Box  c)  depicts  the  HOMO-1  for  AlCp  and  overlap  between  the  A1  p 
and  the  Cp  e\. 


Figure  13.  AlCp  Orbitals 


For  the  simplest  system,  AlCp,  we  see  results  that  provide  good  agreement 
between  the  DFTB  molecular  orbitals  and  the  corresponding  calculated  molecular 
orbitals  as  calculated  from  DFT  in  previous  research.  The  DFT  calculations  for  the 
molecular  orbitals  were  performed  with  the  B3LYP  DFT  functional  using  the  6-31(g,p) 
basis  set  [12].  The  molecular  orbitals  depicted  here  are  three  of  the  bonding  orbitals  for 
AlCp,  the  HOMO,  HOMO-1,  and  HOMO-7  (Figure  13)  and  correspond  to  those 
calculated  in  [12]. 
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HOMO-1 


HOMO-7 


Figure  14.  AlCp  Energy  Levels 

With  the  energy  levels,  we  observe  in  Figure  14  that  DFTB  underestimates  the 
energy  levels  of  the  orbitals  relative  to  our  DFT  results.  The  degenerate  pairing  of 
molecular  orbital  energies  is  consistent  between  both  DFTB  and  DFT.  The  corresponding 
orbital  energies  are  shifted  down  by  approximately  between  1.0  eV  to  1.3  eV  from  DFT 
to  DFTB. 
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a)  AlCp 


[HOMO 

Box  a)  depicts  the  HOMO  for  AlCp*  and  its  nonbonding  interaction  between  the  Cp  ci\ 
and  the  A1  sp  as  well  as  an  additional  contribution  from  the  added  methyl  groups.  Box  b) 
depicts  the  HOMO-5  for  AlCp*  and  its  Cp  (t\  bonding  interaction  with  the  A1  sp  as  well 
as  the  additional  contribution  from  the  added  methyl  groups. 

Figure  15.  AlCp*  Orbitals 

For  the  AlCp*  system  (Figure  15),  we  can  see  good  parallels  with  the  AlCp 
system’s  molecular  orbitals.  Here  the  HOMO  is  generally  in  good  agreement  with  the 
HOMO  for  the  AlCp  system,  though  with  notable  differences  accounting  for  the  addition 
of  the  methyl  groups  to  the  outside  of  the  Cp  ring.  The  HOMO-5  is  analogous  to  the 
HOMO-7  for  AlCp,  again  with  additional  orbital  contribution  corresponding  to  the  added 
methyl  groups.  However,  no  analogous  molecular  orbital  for  the  HOMO-1  in  AlCp  could 
be  found  for  AlCp*. 
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Figure  16.  AlCp*  Energy  Levels 

Once  again,  our  results  (Figure  16)  indicate  that  DFTB  underestimates  the  energy 
levels  of  the  orbitals  relative  to  our  DFT  results  while  the  degenerate  pairing  remains 
consistent.  The  corresponding  orbital  energies  are  shifted  down  by  approximately 
between  0.5  eV  to  0.7  eV  from  DFT  to  DFTB. 
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Box  a)  depicts  the  HOMO  for  AlCp*Pl  and  overlap  between  the  A1  p  and  the  Cp  e\  as 
well  as  an  additional  contribution  from  the  attached  Propyl  group.  Box  b)  depicts  the 
HOMO-2  for  AlCp*Pl  its  nonbonding  interaction  between  the  Cp  ii\  and  the  A1  sp  as  well 
as  the  additional  contribution  from  the  attached  Propyl  group.  Note  that  the  relative 
ordering  of  the  HOMO  and  HOMO-2  for  this  system  is  in  reverse  order  from  the 
analogous  systems  in  AlCp. 

Figure  17.  AlCp*Pl  Orbitals 


For  the  AlCp*Pr  system  (Figure  17),  we  again  see  orbitals  analogous  to  those 
found  in  AlCp.  However,  there  is  no  system  analogous  to  the  HOMO-7  of  AlCp. 
Additionally,  the  HOMO  of  AlCp*Pr  is  analogous  to  the  HOMO-1  of  AlCp,  and  the 
HOMO-2  of  AlCp*Pr  is  analogous  to  the  HOMO  of  AlCp,  indicating  that  the  relative 
ordering  of  those  bonding  orbitals  is  switched.  A  possible  explanation  for  this  behavior 
appears  in  Figure  18. 
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Figure  18.  AlCp*Pr  Energy  Levels 

Figure  18  shows  a  possible  explanation  for  the  change  in  ordering  of  the  HOMO 
and  HOMO-2  orbitals  of  AlCp*Pr  relative  to  AlCp’s  analogous  orbitals.  The  energy 
levels  of  the  three  highest  occupied  orbitals  in  our  DFTB  calculation  are  very  close, 
nearly  triply  degenerate.  The  general  trend  for  underestimating  the  energy  levels  and 
degenerate  pairing  continues  for  our  third  monomer.  The  corresponding  orbital  energies 
are  shifted  down  by  approximately  between  0.3  eV  to  0.6  eV  from  DFT  to  DFTB. 
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Figure  19.  Al4Cp4  Orbitals 

The  applicability  of  waveplot  calculated  molecular  orbitals  rapidly  breaks  down 
for  more  complex  systems.  For  the  simplest  tetramer,  Al4Cp4,  we  can  see  in  Figure  19 
that  the  HOMO- 1 1  orbital  here  does  not  match  the  equivalent  bonding  orbital  obtained 
through  DFT.  Most  notably,  there  is  considerable  asymmetry  between  the  lower  3 
ligands,  as  well  as  significant  distortion  around  the  A1  core.  We  believe  that  the  minimal 
basis  set  used  in  DFTB,  while  providing  us  with  dramatic  improvements  in  calculation 
speed,  leads  to  a  significant  loss  of  accuracy  for  more  complex  systems,  at  least  without  a 
more  thorough  re-parametrization. 
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Figure  20.  Al4Cp4  Energy  Levels 
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0  - 

n  - 

-0.012 
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-0.050 

-0.056 

_n  1:70 
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-0.653 

-0.973 
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-1.014 
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-1.223 

-1.334 

-1.353 

-1.528 

-1.541 

1  -l.DUb 

-1.603 

-1.576 

-1.635 

Figure  21.  AUCp *4  Energy  Levels 
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A  general  trend  is  established  for  the  energy  levels  of  the  tetramer  family  of 
systems  as  presented  in  Figure  20,  Figure  21,  and  Figure  22.  The  first  two  energy  levels 
below  the  HOMO  are  in  relatively  good  agreement  between  DFT  and  DFTB.  However, 
the  subsequent  3  energy  levels  are  significantly  underestimated,  with  the  difference  from 
DFT  and  DFTB  ranging  from  approximately  0.6  to  0.8  eV.  The  next  two  energy  levels 
have  a  difference  with  a  much  wider  range,  with  the  underestimation  from  DFT  to  DFTB 
starting  at  approximately  0.4  eV  for  A^Cp*  4  and  going  to  approximately  0.8  eV  for 
AI4CP4.  The  subsequent  three  energy  levels  are  underestimated  by  approximately  0.4  eV 
for  AI4CP4  going  from  DFT  to  DFTB,  and  is  in  generally  good  agreement  between  the 
two  methods  for  Al4Cp*4  and  Al4Cp*Pr4. 
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Figure  23.  Al8Cp*4  Energy  Levels 


The  Al8Cp*4  is  the  first  indication  of  problems  in  matching  the  degeneracy  of 
energy  levels.  As  we  see  in  Figure  23,  after  the  first  two  energy  levels  below  the  HOMO, 
the  groupings  of  energy  levels  breaks  down.  Overall,  we  see  that  the  performance  of  the 
energy  level  calculations  becomes  more  and  more  problematic  as  the  complexity  of  our 
system  increases,  particularly  as  we  introduce  additional  A1  atoms  to  the  system.  As  seen 
with  the  orbital  plots,  we  believe  that  the  minimal  basis  set  requires  us  to  conduct  a  re- 
parametrization  of  some  key  interactions  to  improve  our  performance  for  more  complex 
systems. 

D.  MOLECULAR  DYNAMICS 

Recently  our  group  has  performed  large-scale  ab  initio  molecular  dynamics  of  the 
oxidation  of  ALCp  4  using  a  DFT-based  Car-Parrinello  method.  This  simulation,  which 
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several  months  on  the  DOD  supercomputers,  was  able  to  successfully  predict  many  of  the 
major  decomposition  products  of  this  cluster  and  was  validated  by  recent  experimental 
studies.  In  this  section  we  discuss  efforts  to  use  this  newly  developed  potential  to  more 
efficiently  simulate  chemistry  in  Al/Cp  clusters.  We  performed  molecular  dynamics  in 

>(t 

DFTB+  with  a  canonical  (NVT)  ensemble  for  the  AI4CP4  system  and  an  AI3O2CP  3 
system,  which  is  the  initial  expected  reaction  product  of  the  AI4CP4  cluster  with  an  O2 
molecule.  An  Andersen  thermostat  was  used,  and  the  system  was  advanced  using  a 
velocity  Verlet  algorithm  with  a  timestep  of  0.05  femtoseconds. 

The  purpose  of  the  AI4CP4  simulation  was  to  observe  an  equilibration  of  the 
system  at  a  temperature  of  300K,  and  determine  if  the  system  is  stable  against  distortions 
when  “annealed”  with  a  simulated  temperature.  In  Figure  24a  below,  we  see  the  initial 
state  of  the  system.  From  this  starting  configuration,  we  see  that  at  5  ps,  in  Figure  24b, 
one  of  the  ligands  has  exhibited  significant  movement  from  its  starting  position  in  the 
tetramer.  By  the  termination  of  the  mn  at  25  ps,  we  see  that  another  ligand,  near  the 
bottom  left  of  the  view  in  Figure  24c,  demonstrates  movement  from  its  starting  position 
as  well.  This  suggests  that  the  improved  structures  with  the  new  potential  may  be  in  a 
local  minimum,  and  that  if  given  sufficient  thermal  energy  we  may  again  find  distortion 
away  from  the  known  Al-Cp  bond  configuration.  In  future  work  a  simulated  annealing 
run  using  constant  temperature  MD  is  recommended  as  an  improved  check  on  the 
repulsive  potential  fits. 
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a)  AI4Cp4 

b) 

c) 

ft 

X 

\ 

t  =  0  ps 

t  =  5  ps 

t  =  25  ps 

Box  a)  depicts  the  initial  state  of  the  system,  Box  b)  depicts  the  system  at  5  ps  into  the 
run,  and  Box  c)  depicts  the  system  at  the  termination  of  the  run  at  25ps.  Note  the  shifting 
in  the  top  ligand  position  from  Box  a)  to  b),  and  then  the  shifting  of  the  bottom  left  ligand 
position  from  Box  b)  to  c). 

Figure  24.  AI4CP4  Molecular  Dynamics 

* 

We  also  simulated  the  molecular  dynamics  of  the  AI3O2CP  3  reacted  product.  In 
Figure  25a  we  can  see  the  starting  configuration  of  the  system.  Figure  25b  depicts  the 
system  after  5ps  of  simulation  time,  and  we  see  that  the  ligands  have  begun  to  shift 
slightly.  In  Figure  25c,  after  10.7ps,  we  see  an  encouraging  result  in  which  a  hydride 
transfer  from  one  of  the  Cp  methyl  groups  to  an  A1  center.  This  reaction  also  occurs  in 
the  DFT  simulations,  and  there  is  indirect  evidence  for  it  in  recent  hot-write  T-jump/mass 
spec  experiments  performed  at  University  of  Maryland  on  small  quantities  of  the 
tetramer.  In  Figure  25d,  after  15ps,  we  can  see  a  second  hydride  has  broken  off  to  bond 
with  another  A1  atom.  In  Figure  25e,  a  third  hydride  transfer  has  occurred.  In  Figure  25f, 
the  ligands  appear  to  have  detached  from  the  core,  and  the  core  has  straightened  into  a 
new  geometry.  Though  these  results  do  not  exactly  reflect  the  DFT  results,  this  DFTB 
simulation  provides  a  very  reasonable  approximation  of  the  reaction  geometry  seen  in 
DFT,  and  captures  some  of  the  difficult  effects  such  as  the  hydride  transfer. 
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a)  AI302Cp3 

t  =  0  ps 

b) 

t  =  5  ps 

c) 

-1  Wli. 

t  =  10.7  ps 

d)  AI302Cp3 

t  =  15  ps 

e) 

t  =  17.25  ps 

f) 

t  =  25  ps 

Box  a)  depicts  the  initial  state  of  the  system.  Box  b)  depicts  the  system  at  5  ps  into  the 
run.  Box  c)  depicts  the  system  at  10.7  ps.  Box  d)  depicts  the  system  at  15ps,  Box  e) 
depicts  the  system  at  17.25  ps,  and  Box  f)  depicts  the  system  at  the  termination  of  the  run 
at  25ps.  Note  the  H  atom  transfer  in  Box  c),  another  H  atom  transfer  in  Box  d),  and  a 
third  H  atom  transfer  in  Box  e).  Box  f)  depicts  the  final  configuration,  where  the  ligands 
appear  to  have  detached  from  the  core,  which  has  reconfigured  into  a  new  geometry. 

Figure  25.  AI3O2CP3  Molecular  Dynamics 
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V.  CONCLUSIONS 


A.  PERFORMANCE 

We  have  demonstrated  that  it  is  possible  to  significantly  improve  the  performance 
of  DFTB  for  the  aluminum  cyclopentadienyl  family  of  clusters  with  the  re- 
parametrization  of  just  a  single  element  pair  interaction.  Though  our  new  parametrization 
demonstrates  some  shortfalls  in  more  complex  systems,  we  were  looking  to  gain  the  most 
return  on  effort  by  targeting  the  most  egregious  shortcomings,  and  that  this  is  far  from  a 
full  and  thorough  re-parametrization. 

The  geometry  optimization  was  the  most  significant  improvement  provided  by  the 
Al-C  re-parametrization.  In  particular,  we  were  able  to  correct  the  gross  geometry 
distortions  produced  by  the  stock  matsci-0-3  parameter  set  in  the  monomers.  With  a 
testing  process  to  fine  tune  the  re-parametrization,  we  were  able  to  extend  this 
improvement  to  more  complex  geometries,  such  as  the  tetramers.  While  we  saw  a  loss  in 
accuracy  as  we  moved  to  even  more  complex  geometries,  the  distortions  did  not  appear 
to  progressively  grow  worse  as  the  system  complexity  grew  and  the  optimized  geometries 
were  still  reasonably  accurate. 

Molecular  orbitals  results  were  less  ideal.  Though  we  saw  good  accuracy  in  the 
orbitals  for  monomers,  increasing  complexity  to  tetramers  led  to  grossly  inaccurate 
distortions.  In  addition,  energy  levels  were  not  in  good  agreement  with  DFT  results. 
Molecular  dynamics  showed  considerable  promise  for  future  research.  Though  not  in 
perfect  agreement  with  DFT,  the  reaction  process  for  the  reacted  product  simulated  with 
reparametrized  DFTB  captures  some  of  the  key  features  seen  in  a  DFT  result. 

B.  FUTURE  WORK 

The  tremendous  reduction  in  the  computational  requirement  provided  by  DFTB 

for  these  systems  is  a  powerful  incentive  for  further  refinement  of  this  method.  The 

geometry  optimization  of  the  AlsoCp*i2  system  was  performed  in  under  24  hours  with 

DFTB  on  a  single  processor,  a  process  that  can  take  thousands  of  CPU  hours  with  higher 

level  levels  of  DFT.  Similar  gains  in  processing  time  were  seen  in  the  molecular 
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dynamics  simulations,  with  each  run  taking  under  72  hours,  in  contrast  to  DFT  run  times 
measured  in  months. 

For  future  work,  it  is  recommended  that  a  more  thorough  re-parametrization  be 
conducted  based  on  performance  metrics  beyond  geometry  optimization.  There  are 
several  other  parameters  that  can  be  adjusted  in  a  subsequent  attempt  to  improve 
performance.  To  improve  the  geometry  results,  additional  improvements  to  accuracy 
could  be  gained  by  reparametrizing  other  interactions  in  the  systems.  Additionally, 
refinements  to  the  band  structure  energy  or  Coulomb  energy  through  adjusting  the  DFT 
method  used  or  parameters  such  as  the  Hubbard  parameter  or  electronegativity  could 
improve  the  electronic  structure  calculations  and  reduce  the  magnitude  of  “correction” 
characterized  in  the  repulsive  potential.  Lastly,  the  Al-C  interaction  demonstrated  that 
geometry  results  were  highly  sensitive  to  even  small  changes  in  the  parametrization. 
Therefore,  a  spline  fit  may  offer  better  results  by  allowing  adjustments  to  more  easily 
focus  upon  the  interaction  region  of  interest. 
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APPENDIX  A.  FUNDAMENTAL  CONSTANTS 


1  Ha  =  27.21 14  eV 
1  Bohr  =  0.5292  A 
lu  =  1.6605  x  10  27  kg 
Time:  1.0327  fs 

With  Hartree  atomic  units,  the  following  constants  are  unity  by  definition: 
me  (electron  mass) 
e  (elementary  charge) 
h  (reduced  Planck’s  constant) 
ke  (Coulomb’s  constant) 
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APPENDIX  B.  ALUMINUM-CARBON  REPARAMETRIZATION 


0.00,  0.0352  -0.067  0.0585  -0.0273  0.0069  -0.0010  0.0001  0,  6.81,  10*0.0 

The  above  text  is  line  2  of  the  Al-C  .skf  file  which  contains  the  re-parametrization 
of  the  Al-C  interactions.  This  line  is  manipulated  identically  in  the  similar  C-Al  .skf  file. 
The  formatting  used  for  this  interaction  is  the  hetero-nuclear  case  in  the  DFTB  skf 
documentation. 

The  first  portion  (0.00,)  is  a  placeholder,  and  is  used  for  atomic  mass  in  homo- 
nuclear  cases.  The  second  portion  (0.0352  -0.067  0.0585  -0.0273  0.0069  -0.0010  0.0001 
0,)  is  the  polynomial,  with  the  input  values  reflecting  the  cn  coefficient  values.  The  next 
portion  (6.81,)  is  the  cutoff  radius,  and  these  values  are  specified  in  Bohr.  The  final 
portion  (10*0.0)  is  just  a  placeholder. 

This  line  of  the  file  would  be  modified  to  adjust  the  repulsive  potential  portion  of 
the  DFTB  parametrization. 
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